"""
Phase 7: Five Extension Tests
===============================
1. OECD bundle decomposition — pension spending, health expenditure as
   specific mechanisms behind the OECD savings attenuation
2. Pooled monetary union — EMU+CFA+ECCU+CMA as single union dummy for power
3. Time stability — does the coefficient surface shift pre/post GFC?
4. Out-of-sample validation — estimate on pre-2010, predict post-2010
5. KAOPEN channel test — does openness matter for channel (trade vs income)
   rather than magnitude?

At the end: honest assessment of what's real vs what's noise.
"""

import pandas as pd
import numpy as np
from pathlib import Path
import sys
import warnings
warnings.filterwarnings('ignore')

PROJECT_DIR = Path(__file__).resolve().parent.parent
ROOT_DIR = PROJECT_DIR.parent
sys.path.insert(0, str(ROOT_DIR / "multilateral" / "src"))
from model import PanelGLS

DATA = PROJECT_DIR / "data" / "processed"
OUT_TABLES = PROJECT_DIR / "output" / "tables"

Z_VARS = ['Z_1', 'Z_2', 'Z_3']


def stars(p):
    if p < 0.01: return '***'
    if p < 0.05: return '**'
    if p < 0.1: return '*'
    return ''


def fmt(val, p):
    return f"{val:.3f}{stars(p)}"


def run_gls(df, y_var, x_vars):
    cols = [y_var] + x_vars + ['iso3', 'year']
    sub = df[cols].dropna()
    if len(sub) < 50:
        return None
    gls = PanelGLS()
    try:
        gls.fit(sub[y_var].values, sub[x_vars].values,
                sub['iso3'].values, sub['year'].values)
    except Exception:
        return None
    result = {'n_obs': gls.n_obs, 'n_countries': gls.n_countries,
              'r_squared': gls.r_squared}
    for i, var in enumerate(x_vars):
        result[f'{var}_coef'] = gls.beta[i]
        result[f'{var}_se'] = gls.se[i]
        result[f'{var}_p'] = gls.pvalues[i]
    return result


# ── Union definitions ─────────────────────────────────────────────────

EUROZONE_JOIN = {
    'AUT': 1999, 'BEL': 1999, 'FIN': 1999, 'FRA': 1999, 'DEU': 1999,
    'IRL': 1999, 'ITA': 1999, 'LUX': 1999, 'NLD': 1999, 'PRT': 1999,
    'ESP': 1999, 'GRC': 2001, 'SVN': 2007, 'CYP': 2008, 'MLT': 2008,
    'SVK': 2009, 'EST': 2011, 'LVA': 2014, 'LTU': 2015,
}

CFA_JOIN = {
    'BEN': 1960, 'BFA': 1960, 'CIV': 1960, 'MLI': 1984,
    'NER': 1960, 'SEN': 1960, 'TGO': 1960, 'GNB': 1997,
    'CMR': 1960, 'CAF': 1960, 'TCD': 1960, 'COG': 1960,
    'GNQ': 1985, 'GAB': 1960,
}

ECCU_JOIN = {
    'ATG': 1976, 'DMA': 1976, 'GRD': 1976, 'KNA': 1976,
    'LCA': 1976, 'VCT': 1976, 'MSR': 1976, 'AIA': 1976,
}

CMA_JOIN = {
    'ZAF': 1974, 'LSO': 1974, 'SWZ': 1974, 'NAM': 1990,
}


def add_union_flags(df):
    """Add monetary union indicators."""
    df['in_emu'] = 0
    df['in_cfa'] = 0
    df['in_eccu'] = 0
    df['in_cma'] = 0
    df['in_any_union'] = 0

    for iso, yr in EUROZONE_JOIN.items():
        df.loc[(df['iso3'] == iso) & (df['year'] >= yr), 'in_emu'] = 1
    for iso, yr in CFA_JOIN.items():
        df.loc[(df['iso3'] == iso) & (df['year'] >= yr), 'in_cfa'] = 1
    for iso, yr in ECCU_JOIN.items():
        df.loc[(df['iso3'] == iso) & (df['year'] >= yr), 'in_eccu'] = 1
    for iso, yr in CMA_JOIN.items():
        df.loc[(df['iso3'] == iso) & (df['year'] >= yr), 'in_cma'] = 1

    df['in_any_union'] = ((df['in_emu'] == 1) | (df['in_cfa'] == 1) |
                          (df['in_eccu'] == 1) | (df['in_cma'] == 1)).astype(int)
    return df


def main():
    print("Phase 7: Extension Tests")
    print("=" * 70)

    df = pd.read_csv(DATA / "unified_panel.csv")
    df = add_union_flags(df)
    print(f"Panel: {len(df)} obs, {df['iso3'].nunique()} countries")
    print(f"Union obs: EMU={df['in_emu'].sum()}, CFA={df['in_cfa'].sum()}, "
          f"ECCU={df['in_eccu'].sum()}, CMA={df['in_cma'].sum()}, "
          f"Any={df['in_any_union'].sum()}")

    # Create union interaction terms
    for z in Z_VARS:
        df[f'{z}_x_union'] = df[z] * df['in_any_union']

    # ══════════════════════════════════════════════════════════════════
    # TEST 1: OECD Bundle Decomposition
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("TEST 1: OECD Bundle Decomposition")
    print("What specifically attenuates the OECD savings channel?")
    print("-" * 50)

    # Create pension/health interaction terms
    # Pension spending as % GDP — use as continuous moderator
    df['has_pension'] = df['pension_spending_gdp'].notna().astype(int)
    # Standardize pension spending for interaction
    pen_mean = df.loc[df['pension_spending_gdp'].notna(), 'pension_spending_gdp'].mean()
    pen_std = df.loc[df['pension_spending_gdp'].notna(), 'pension_spending_gdp'].std()
    df['pension_std'] = (df['pension_spending_gdp'] - pen_mean) / pen_std
    df['pension_std'] = df['pension_std'].fillna(0)
    df['pension_high'] = (df['pension_spending_gdp'] > df['pension_spending_gdp'].median()).astype(int)
    df['pension_high'] = df['pension_high'].fillna(0)

    # Health spending
    he_mean = df.loc[df['health_exp_gdp'].notna(), 'health_exp_gdp'].mean()
    he_std = df.loc[df['health_exp_gdp'].notna(), 'health_exp_gdp'].std()
    df['health_std'] = (df['health_exp_gdp'] - he_mean) / he_std
    df['health_std'] = df['health_std'].fillna(0)

    for z in Z_VARS:
        df[f'{z}_x_pension'] = df[z] * df['pension_std']
        df[f'{z}_x_penhigh'] = df[z] * df['pension_high']
        df[f'{z}_x_health'] = df[z] * df['health_std']

    test1_results = []
    for dv in ['gross_savings_gdp', 'ca_gdp']:
        print(f"\n  DV: {dv}")

        # M1: Z + Z×OECD (baseline from Phase 2)
        x1 = Z_VARS + ['Z_1_x_oecd', 'Z_2_x_oecd', 'Z_3_x_oecd']
        m1 = run_gls(df, dv, x1)

        # M2: Z + Z×pension_high
        x2 = Z_VARS + ['Z_1_x_penhigh', 'Z_2_x_penhigh', 'Z_3_x_penhigh']
        m2 = run_gls(df, dv, x2)

        # M3: Z + Z×pension_high + Z×OECD
        x3 = Z_VARS + ['Z_1_x_penhigh', 'Z_2_x_penhigh', 'Z_3_x_penhigh',
                        'Z_1_x_oecd', 'Z_2_x_oecd', 'Z_3_x_oecd']
        m3 = run_gls(df, dv, x3)

        # M4: Z + Z×pension_std (continuous)
        x4 = Z_VARS + ['Z_1_x_pension', 'Z_2_x_pension', 'Z_3_x_pension']
        m4 = run_gls(df, dv, x4)

        # M5: Z + Z×health_std
        x5 = Z_VARS + ['Z_1_x_health', 'Z_2_x_health', 'Z_3_x_health']
        m5 = run_gls(df, dv, x5)

        # M6: Z + Z×pension + Z×health + Z×OECD (full horse race)
        x6 = Z_VARS + ['Z_1_x_pension', 'Z_2_x_pension', 'Z_3_x_pension',
                        'Z_1_x_health', 'Z_2_x_health', 'Z_3_x_health',
                        'Z_1_x_oecd', 'Z_2_x_oecd', 'Z_3_x_oecd']
        m6 = run_gls(df, dv, x6)

        models = [('M1: +OECD', m1, 'Z_1_x_oecd'),
                  ('M2: +Pension high', m2, 'Z_1_x_penhigh'),
                  ('M3: +Pension+OECD', m3, 'Z_1_x_penhigh'),
                  ('M4: +Pension(cont)', m4, 'Z_1_x_pension'),
                  ('M5: +Health', m5, 'Z_1_x_health'),
                  ('M6: Full', m6, 'Z_1_x_pension')]

        for name, m, key in models:
            if m is None:
                print(f"    {name}: failed")
                continue
            z1 = fmt(m['Z_1_coef'], m['Z_1_p'])
            ki = fmt(m[f'{key}_coef'], m[f'{key}_p'])
            print(f"    {name}: Z₁={z1}, {key}={ki}, N={m['n_obs']}, R²={m['r_squared']:.3f}")
            row = {'dv': dv, 'model': name, 'Z1': m['Z_1_coef'], 'Z1_p': m['Z_1_p'],
                   'n': m['n_obs'], 'r2': m['r_squared']}
            # Store the key interaction
            row['key_var'] = key
            row['key_coef'] = m[f'{key}_coef']
            row['key_p'] = m[f'{key}_p']
            # Store OECD if present
            if f'Z_1_x_oecd_coef' in m:
                row['oecd_coef'] = m['Z_1_x_oecd_coef']
                row['oecd_p'] = m['Z_1_x_oecd_p']
            test1_results.append(row)

    # ══════════════════════════════════════════════════════════════════
    # TEST 2: Pooled Monetary Union
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("TEST 2: Pooled Monetary Union Interaction")
    print(f"Any-union obs: {df['in_any_union'].sum()} "
          f"({df.loc[df['in_any_union']==1, 'iso3'].nunique()} countries)")
    print("-" * 50)

    test2_results = []
    for dv in ['ca_gdp', 'gross_savings_gdp', 'gross_investment_gdp', 'nfa_gdp']:
        if dv not in df.columns:
            continue

        # M1: Z + Z×EMU only (from Phase 2)
        x1 = Z_VARS + ['Z_1_x_emu', 'Z_2_x_emu', 'Z_3_x_emu']
        m1 = run_gls(df, dv, x1)

        # M2: Z + Z×any_union (pooled)
        x2 = Z_VARS + ['Z_1_x_union', 'Z_2_x_union', 'Z_3_x_union']
        m2 = run_gls(df, dv, x2)

        print(f"\n  {dv}:")
        if m1:
            print(f"    EMU only:   Z₁×EMU={fmt(m1['Z_1_x_emu_coef'], m1['Z_1_x_emu_p'])}, "
                  f"N={m1['n_obs']}")
        if m2:
            print(f"    Any union:  Z₁×union={fmt(m2['Z_1_x_union_coef'], m2['Z_1_x_union_p'])}, "
                  f"N={m2['n_obs']}")
            test2_results.append({
                'dv': dv,
                'emu_coef': m1['Z_1_x_emu_coef'] if m1 else np.nan,
                'emu_p': m1['Z_1_x_emu_p'] if m1 else np.nan,
                'union_coef': m2['Z_1_x_union_coef'],
                'union_p': m2['Z_1_x_union_p'],
                'n_union': m2['n_obs'],
            })

    # ══════════════════════════════════════════════════════════════════
    # TEST 3: Time Stability (Pre/Post GFC)
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("TEST 3: Time Stability of Coefficient Surface")
    print("-" * 50)

    pre_gfc = df[df['year'] <= 2007]
    post_gfc = df[df['year'] >= 2009]
    print(f"  Pre-GFC: {len(pre_gfc)} obs, Post-GFC: {len(post_gfc)} obs")

    test3_results = []
    for dv in ['ca_gdp', 'gross_savings_gdp']:
        if dv not in df.columns:
            continue

        # Baseline and key interactions for each period
        for period_name, period_df in [('Full', df), ('Pre-GFC', pre_gfc), ('Post-GFC', post_gfc)]:
            # Baseline
            m0 = run_gls(period_df, dv, Z_VARS)
            # Income interaction
            x_inc = Z_VARS + ['Z_1_x_high', 'Z_2_x_high', 'Z_3_x_high',
                              'Z_1_x_low', 'Z_2_x_low', 'Z_3_x_low']
            m_inc = run_gls(period_df, dv, x_inc)
            # OECD interaction
            x_oecd = Z_VARS + ['Z_1_x_oecd', 'Z_2_x_oecd', 'Z_3_x_oecd']
            m_oecd = run_gls(period_df, dv, x_oecd)

            if m0 and m_inc and m_oecd:
                print(f"\n  {dv} [{period_name}]:")
                print(f"    Baseline Z₁={fmt(m0['Z_1_coef'], m0['Z_1_p'])}, N={m0['n_obs']}")
                print(f"    Z₁×high={fmt(m_inc['Z_1_x_high_coef'], m_inc['Z_1_x_high_p'])}, "
                      f"Z₁×low={fmt(m_inc['Z_1_x_low_coef'], m_inc['Z_1_x_low_p'])}")
                print(f"    Z₁×OECD={fmt(m_oecd['Z_1_x_oecd_coef'], m_oecd['Z_1_x_oecd_p'])}")

                test3_results.append({
                    'dv': dv, 'period': period_name,
                    'Z1_base': m0['Z_1_coef'], 'Z1_base_p': m0['Z_1_p'],
                    'Z1xhigh': m_inc['Z_1_x_high_coef'], 'Z1xhigh_p': m_inc['Z_1_x_high_p'],
                    'Z1xlow': m_inc['Z_1_x_low_coef'], 'Z1xlow_p': m_inc['Z_1_x_low_p'],
                    'Z1xoecd': m_oecd['Z_1_x_oecd_coef'], 'Z1xoecd_p': m_oecd['Z_1_x_oecd_p'],
                    'n': m0['n_obs'],
                })

    # ══════════════════════════════════════════════════════════════════
    # TEST 4: Out-of-Sample Validation
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("TEST 4: Out-of-Sample Validation")
    print("Estimate surface on pre-2010, predict post-2010 subsample effects")
    print("-" * 50)

    train = df[df['year'] <= 2009]
    test = df[df['year'] >= 2010]

    test4_results = []
    for dv in ['ca_gdp', 'gross_savings_gdp']:
        if dv not in df.columns:
            continue

        # Estimate full interaction model on training data
        x_full = Z_VARS + ['Z_1_x_high', 'Z_2_x_high', 'Z_3_x_high',
                           'Z_1_x_low', 'Z_2_x_low', 'Z_3_x_low',
                           'Z_1_x_oecd', 'Z_2_x_oecd', 'Z_3_x_oecd']
        m_train = run_gls(train, dv, x_full)
        if m_train is None:
            continue

        # Predict: implied Z₁ for OECD vs non-OECD, high vs low income
        z1_base_train = m_train['Z_1_coef']
        z1_high_train = m_train['Z_1_coef'] + m_train['Z_1_x_high_coef']
        z1_low_train = m_train['Z_1_coef'] + m_train['Z_1_x_low_coef']
        z1_oecd_train = m_train['Z_1_coef'] + m_train['Z_1_x_high_coef'] + m_train['Z_1_x_oecd_coef']

        # Actual: estimate Z₁ on test-period subsamples
        subsamples = {
            'Full': test,
            'High income': test[test['income_high'] == 1],
            'Low income': test[test['income_low'] == 1],
            'OECD': test[test['is_oecd'] == 1],
        }

        predicted = {
            'Full': z1_base_train,
            'High income': z1_high_train,
            'Low income': z1_low_train,
            'OECD': z1_oecd_train,
        }

        print(f"\n  {dv}:")
        print(f"  {'Subsample':20s} {'Predicted':>10s} {'Actual':>12s} {'Match?':>8s}")
        for sname, sdf in subsamples.items():
            m_actual = run_gls(sdf, dv, Z_VARS)
            if m_actual is None:
                continue
            pred = predicted[sname]
            actual = m_actual['Z_1_coef']
            actual_p = m_actual['Z_1_p']
            # Match = same sign and within 2x magnitude
            same_sign = np.sign(pred) == np.sign(actual) if pred != 0 and actual != 0 else abs(pred - actual) < 10
            ratio = abs(actual / pred) if pred != 0 else float('inf')
            match = same_sign and 0.25 < ratio < 4.0
            print(f"  {sname:20s} {pred:10.1f} {fmt(actual, actual_p):>12s} {'YES' if match else 'NO':>8s}")
            test4_results.append({
                'dv': dv, 'subsample': sname,
                'predicted': pred, 'actual': actual, 'actual_p': actual_p,
                'match': match, 'n': m_actual['n_obs'],
            })

    # ══════════════════════════════════════════════════════════════════
    # TEST 5: KAOPEN Channel Test
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("TEST 5: KAOPEN — Channel Routing vs Magnitude")
    print("Does openness change WHICH channel demographics use?")
    print("-" * 50)

    # Need CA components — check what's available
    ca_components = ['trade_openness']  # proxy; actual components would need WDI
    # Use savings and investment as the S-I decomposition
    test5_results = []
    for channel_dv, label in [('gross_savings_gdp', 'Savings'),
                               ('gross_investment_gdp', 'Investment'),
                               ('ca_gdp', 'CA')]:
        if channel_dv not in df.columns:
            continue

        # Z + Z×KAOPEN (continuous)
        x = Z_VARS + ['Z_1_x_kaopen', 'Z_2_x_kaopen', 'Z_3_x_kaopen']
        m = run_gls(df, channel_dv, x)
        if m:
            print(f"  {label}: Z₁={fmt(m['Z_1_coef'], m['Z_1_p'])}, "
                  f"Z₁×KAOPEN={fmt(m['Z_1_x_kaopen_coef'], m['Z_1_x_kaopen_p'])}, "
                  f"N={m['n_obs']}")
            test5_results.append({
                'channel': label,
                'Z1_base': m['Z_1_coef'], 'Z1_base_p': m['Z_1_p'],
                'Z1xKA': m['Z_1_x_kaopen_coef'], 'Z1xKA_p': m['Z_1_x_kaopen_p'],
                'n': m['n_obs'],
            })

        # Also test KAOPEN saturated
        x2 = Z_VARS + ['Z_1_x_ksat', 'Z_2_x_ksat', 'Z_3_x_ksat']
        m2 = run_gls(df, channel_dv, x2)
        if m2:
            print(f"  {label} (sat): Z₁={fmt(m2['Z_1_coef'], m2['Z_1_p'])}, "
                  f"Z₁×KAsat={fmt(m2['Z_1_x_ksat_coef'], m2['Z_1_x_ksat_p'])}")

    # ══════════════════════════════════════════════════════════════════
    # WRITE ALL OUTPUT TABLES
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("Writing output tables...")

    # Test 1
    t1_df = pd.DataFrame(test1_results)
    with open(OUT_TABLES / "phase7_test1_oecd_bundle.md", 'w') as f:
        f.write("# Test 1: OECD Bundle Decomposition\n\n")
        f.write("What institutional feature attenuates the OECD savings channel?\n\n")
        for dv in ['gross_savings_gdp', 'ca_gdp']:
            sub = t1_df[t1_df['dv'] == dv]
            if sub.empty:
                continue
            f.write(f"\n## {dv}\n\n")
            f.write("| Model | Z₁ | Key Interaction | OECD (if present) | N | R² |\n")
            f.write("|---|---|---|---|---|---|\n")
            for _, row in sub.iterrows():
                z1 = fmt(row['Z1'], row['Z1_p'])
                ki = fmt(row['key_coef'], row['key_p'])
                oecd = fmt(row['oecd_coef'], row['oecd_p']) if pd.notna(row.get('oecd_coef')) else '--'
                f.write(f"| {row['model']} | {z1} | {ki} | {oecd} | "
                        f"{int(row['n'])} | {row['r2']:.3f} |\n")
        f.write("\n*Pension spending standardized (mean=0, SD=1). "
                "Pension high = above median pension/GDP.*\n")
    print("  Wrote: phase7_test1_oecd_bundle.md")

    # Test 2
    with open(OUT_TABLES / "phase7_test2_pooled_union.md", 'w') as f:
        f.write("# Test 2: Pooled Monetary Union Interaction\n\n")
        f.write(f"Pooled union = EMU ({sum(1 for v in EUROZONE_JOIN)}) + "
                f"CFA ({sum(1 for v in CFA_JOIN)}) + "
                f"ECCU ({sum(1 for v in ECCU_JOIN)}) + "
                f"CMA ({sum(1 for v in CMA_JOIN)}) = "
                f"{df['in_any_union'].sum()} obs\n\n")
        f.write("| DV | Z₁×EMU only | p | Z₁×Any Union | p | N |\n")
        f.write("|---|---|---|---|---|---|\n")
        for row in test2_results:
            emu = f"{row['emu_coef']:.3f}" if not pd.isna(row['emu_coef']) else '--'
            emu_p = f"{row['emu_p']:.3f}" if not pd.isna(row['emu_p']) else '--'
            f.write(f"| {row['dv']} | {emu} | {emu_p} | "
                    f"{row['union_coef']:.3f} | {row['union_p']:.3f} | "
                    f"{row['n_union']} |\n")
        f.write("\n*Does pooling all monetary unions provide enough power "
                "to detect the union interaction?*\n")
    print("  Wrote: phase7_test2_pooled_union.md")

    # Test 3
    with open(OUT_TABLES / "phase7_test3_time_stability.md", 'w') as f:
        f.write("# Test 3: Time Stability of the Coefficient Surface\n\n")
        f.write("| DV | Period | Z₁ (base) | Z₁×High | Z₁×Low | Z₁×OECD | N |\n")
        f.write("|---|---|---|---|---|---|---|\n")
        for row in test3_results:
            f.write(f"| {row['dv']} | {row['period']} | "
                    f"{fmt(row['Z1_base'], row['Z1_base_p'])} | "
                    f"{fmt(row['Z1xhigh'], row['Z1xhigh_p'])} | "
                    f"{fmt(row['Z1xlow'], row['Z1xlow_p'])} | "
                    f"{fmt(row['Z1xoecd'], row['Z1xoecd_p'])} | "
                    f"{int(row['n'])} |\n")
        f.write("\n*Pre-GFC: <=2007. Post-GFC: >=2009. "
                "Stable surface = similar interaction signs and magnitudes.*\n")
    print("  Wrote: phase7_test3_time_stability.md")

    # Test 4
    with open(OUT_TABLES / "phase7_test4_oos_validation.md", 'w') as f:
        f.write("# Test 4: Out-of-Sample Validation\n\n")
        f.write("Surface estimated on pre-2010 data; predictions tested against post-2010 actuals.\n\n")
        f.write("| DV | Subsample | Predicted Z₁ | Actual Z₁ | Match? | N |\n")
        f.write("|---|---|---|---|---|---|\n")
        for row in test4_results:
            f.write(f"| {row['dv']} | {row['subsample']} | "
                    f"{row['predicted']:.1f} | {fmt(row['actual'], row['actual_p'])} | "
                    f"{'YES' if row['match'] else 'NO'} | {int(row['n'])} |\n")
        n_match = sum(1 for r in test4_results if r['match'])
        f.write(f"\n*{n_match}/{len(test4_results)} predictions match "
                f"(same sign, within 4x magnitude).*\n")
    print("  Wrote: phase7_test4_oos_validation.md")

    # Test 5
    with open(OUT_TABLES / "phase7_test5_kaopen_channel.md", 'w') as f:
        f.write("# Test 5: KAOPEN — Channel Routing\n\n")
        f.write("Does financial openness change which channel demographics use?\n\n")
        f.write("| Channel | Z₁ (base) | Z₁×KAOPEN | N |\n")
        f.write("|---|---|---|---|\n")
        for row in test5_results:
            f.write(f"| {row['channel']} | {fmt(row['Z1_base'], row['Z1_base_p'])} | "
                    f"{fmt(row['Z1xKA'], row['Z1xKA_p'])} | {int(row['n'])} |\n")
        f.write("\n*KAOPEN is continuous Chinn-Ito index. "
                "Positive Z₁×KAOPEN = openness amplifies demographic channel.*\n")
    print("  Wrote: phase7_test5_kaopen_channel.md")

    # ══════════════════════════════════════════════════════════════════
    # HONEST ASSESSMENT
    # ══════════════════════════════════════════════════════════════════
    print(f"\n{'=' * 70}")
    print("HONEST ASSESSMENT: Signal vs Noise")
    print("=" * 70)

    # Count significant results across all tests
    n_tests_total = 0
    n_sig_05 = 0
    n_sig_10 = 0

    # Test 1: key interactions
    for row in test1_results:
        n_tests_total += 1
        if row['key_p'] < 0.05: n_sig_05 += 1
        if row['key_p'] < 0.10: n_sig_10 += 1

    # Test 2: union interactions
    for row in test2_results:
        n_tests_total += 1
        if row['union_p'] < 0.05: n_sig_05 += 1
        if row['union_p'] < 0.10: n_sig_10 += 1

    # Test 3: time stability (count as confirmatory, not significance tests)
    # Test 4: OOS (match/no match)
    n_oos_match = sum(1 for r in test4_results if r['match'])
    n_oos_total = len(test4_results)

    # Test 5: kaopen channel
    for row in test5_results:
        n_tests_total += 1
        if row['Z1xKA_p'] < 0.05: n_sig_05 += 1
        if row['Z1xKA_p'] < 0.10: n_sig_10 += 1

    print(f"\nNew interaction tests: {n_sig_05}/{n_tests_total} significant at 5%, "
          f"{n_sig_10}/{n_tests_total} at 10%")
    print(f"OOS prediction: {n_oos_match}/{n_oos_total} subsample predictions match")

    # Bonferroni threshold
    bonf_05 = 0.05 / n_tests_total if n_tests_total > 0 else 0.05
    print(f"\nBonferroni threshold (5%): p < {bonf_05:.4f}")
    print(f"Multiple testing concern: running {n_tests_total} interaction tests")
    print(f"Expected false positives at 5%: {n_tests_total * 0.05:.1f}")
    print(f"Actual significant at 5%: {n_sig_05}")
    ratio = n_sig_05 / (n_tests_total * 0.05) if n_tests_total > 0 else 0
    print(f"Ratio (actual/expected): {ratio:.1f}x")
    if ratio > 3:
        print("→ Signal exceeds noise expectation by >3x — likely real")
    elif ratio > 1.5:
        print("→ Modestly above chance — some real signal mixed with noise")
    else:
        print("→ AT OR BELOW chance level — could be entirely spurious")

    # Write assessment
    with open(OUT_TABLES / "phase7_assessment.md", 'w') as f:
        f.write("# Phase 7: Honest Assessment — Signal vs Noise\n\n")
        f.write(f"## Test Statistics\n\n")
        f.write(f"- New interaction tests run: {n_tests_total}\n")
        f.write(f"- Significant at 5%: {n_sig_05}\n")
        f.write(f"- Significant at 10%: {n_sig_10}\n")
        f.write(f"- Expected false positives (5%): {n_tests_total * 0.05:.1f}\n")
        f.write(f"- Ratio (actual/expected): {ratio:.1f}x\n")
        f.write(f"- OOS predictions matching: {n_oos_match}/{n_oos_total}\n\n")
        f.write("## Interpretation\n\n")
        if ratio > 3:
            f.write("Signal substantially exceeds noise expectation. "
                    "The core findings (income moderation, OECD attenuation) "
                    "appear robust to extension.\n")
        elif ratio > 1.5:
            f.write("Signal modestly above chance. Some findings are likely real "
                    "(particularly those that replicate Phase 2-6 patterns), "
                    "but individual marginal results should be treated with caution.\n")
        else:
            f.write("Signal at or below chance level. The extension tests do not "
                    "add reliable new information beyond what Phase 2-6 established. "
                    "Further subdivision risks overfitting.\n")
    print("  Wrote: phase7_assessment.md")

    print(f"\nPhase 7 complete.")


if __name__ == '__main__':
    main()
